G4: Edinaldo de Alencar / Igor Freire / Ramon Araújo / Ricardo Ribeiro
30 de outubro de 2014
dataset_diario é a base com valores de fluxos diários de 2002 a 2009.dataset_mensal é a base com os valores de fluxo mensais de 1992 a 2009.Nota: a informação sobre os dias da semana é redundante, pois pode ser obtida através de:
dia <- weekdays(index(dataset_diario))
Por isso, na leitura da base em dataset_diario, esta coluna foi ignorada.
Conversão dos dados para o formato de série temporal no R.
# Frequency --> número de observações por unidade de tempo
# define a unidade de tempo (e.g. 12: unidade de tempo = ano)
tsMensal <- ts(dataset_mensal, frequency=12, start=1992)
Série temporal diária
tsDiario <- msts(dataset_diario, start=c(2002,1,1), seasonal.periods=c(7, 365.25))
Nota: é importante observar os anos bissextos, 2004 e 2008.
Observar a influência dos feriados. Exemplo, 12 de outubro foi uma quinta-feira em 2006.
A série temporal pode ser decomposta em:
acf(dataset_diario, col="red")
pacf(dataset_diario, col="red")
Realizar prospecções de curto, médio e longo prazo.
Para curto prazo, será utilizada a série temporal com dados diários (tsDiario).
Para médio e longo prazo, será utilizada a série temporal com dados mensais (tsMensal).
Separação das séries temporais em conjuntos de treinamento e de testes.
tsDiario)tsDiarioTrain <- window(tsDiario, end=c(2008,183)) # até 30/06/2008
# De 01/07/2008 a 30/07/2008
tsDiarioTest30Days <- window(tsDiario, start=c(2008,184), end=c(2008,213))
# 45 dias a partir De 1/1/2008
tsDiarioTest45Days <- window(tsDiario, start=c(2008,184), end=c(2008,228))
tsMensal)tsMensalTrain <- window(tsMensal, end=c(2007, 12)) # De Jan 1992 a Dez 2007
tsMensalTest4mth <- window(tsMensal, start=2008, end=c(2008, 4)) # De Jan 2006 a Mar 2006
tsMensalTest6mth <- window(tsMensal, start=2008, end=c(2008, 6)) # De Jan 2006 a Jun 2006
tsMensalTest1yr <- window(tsMensal, start=2008, end=c(2008, 12)) # De Jan 2006 a Jan 2007
tsMensalTest2yr <- window(tsMensal, start=2008, end=c(2009, 12)) # De Jan 2006 a Jan 2008
ver (R. J. Hyndman 2014) e (Leek 2014):
library(forecast)
library(nnet)
Apresenta-se a seguir:
Nota:
Diferentemente do ME, MAE e RMSE, o medidor MAPE é independente da escala.
\[M = \frac{1}{n} \sum \limits_{t=1}^{n} \left| \frac{A_t - F_t}{A_t} \right|,\]
onde \(A_t\) é o valor verdadeiro, \(F_t\) é o valor predito e \(n\) é o número de amostras temporais preditas.
Como desvantagem, o critério MAPE tem comportamento indefinido quando o valor verdadeiro \(A_t=0\) (divisão por zero).
30 dias
## Training set Test set
## Mean 14.656 22.506
## Näive 6.214 4.674
## N-Sazonal 8.720 9.066
## Drift 6.216 4.788
## Regressao 5.868 7.383
## Training set Test set
## Mean 25.58 35.295
## Näive 4.44 5.639
## N-Sazonal 7.68 6.051
## Drift 4.47 6.542
1 ano
## Training set Test set
## Mean 25.582 39.867
## Näive 4.440 5.708
## N-Sazonal 7.680 6.898
## Drift 4.470 4.722
## Regressao 4.741 10.023
## ETS 1.442 2.603
2 anos
## Training set Test set
## Mean 25.582 41.232
## Näive 4.440 6.801
## N-Sazonal 7.680 9.046
## Drift 4.470 5.210
## Regressao 4.741 10.398
## ETS 1.442 2.964
Em um modelo de auto-regressão, prevemos a variável de interesse usando uma combinação linear dos valores passados da variável.
Assim, um modelo auto-regressivo de ordem \(p\), denotado por \(AR(p)\), pode ser escrito como:
\(y_t = c + \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + \dots + \varphi_p y_{t-p} + e_t,\)
onde \(c\) é uma constante, \(e_t\) o ruído branco, e, os valores defasados de \(y_t\) são os preditores.
Para um modelo \(AR(1)\), tem-se \(\varphi_1\in(-1,1)\)
Para um modelo \(AR(2)\), tem-se \(\varphi_1,\varphi_2\in(-1,1)\) com \(\varphi_1 + \varphi_2 < 1\) e \(\varphi_2 - \varphi_1 < 1\)
Ao invés de usar valores passados da variável de previsão em uma regressão, o modelo de médias-móveis usa erros de previsão do passado em um modelo de regressão-like.
\(y_t = c + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \dots + \theta_q e_{t-q}.\)
Esse modelo tem ordem \(q\) e é denotado por \(MA(q)\).
Para um modelo \(MA(1)\), tem-se \(\theta_1\in(-1,1)\)
Para um modelo \(MA(2)\), tem-se \(\theta_1,\theta_2\in(-1,1)\) com \(\theta_1 + \theta_2 > -1\) e \(\theta_1 - \theta_2 < 1\)
É possível escrever o modelo \(AR(p)\) como um modelo \(MA(\infty)\). Em particular para um modelo \(AR(1)\), vem:
\(y_t = \varphi_1 y_{t-1} + e_t\)
\(y_t = \varphi_1 ( \varphi_1 y_{t-2} + e_{t-1} ) + e_t\)
\(y_t = \varphi_1^2 y_{t-2} + \varphi_1 e_{t-1} + e_t\)
\(y_t = \varphi_1^3 y_{t-3} + \varphi_1^2 e_{t-2} + \varphi_1 e_{t-1} + e_t\)
e assim por diante. Onde \(\varphi_1\in(-1,1)\)
ARIMA (Autoregressive Integrated Moving Average)
Método ARMA, acrescido da etapa de integração (I).
Caso a série temporal seja estacionária, o método se reduz ao ARMA.
\[y(t) = \sum_{i=1}^{p}\varphi_ix(t-i) + \sum_{j=1}^{q}\theta_j\epsilon(t-j)\]
\(ARIMA(p,d,q)\), onde:
p = número de valores passados da variável (AR)
d = número de diferenças não-sazonais
q = número de rerros de predição passados (MA)
mesesAPrever <- 4
Modelo:
arima_model_mensal <- auto.arima(tsMensalTrain)
Predição:
fcast_arima_4mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_4mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest4mth, col="red")
Acurácia:
accuracy(fcast_arima_4mth, tsMensalTest4mth)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 716.4 11391 8353 0.115 1.625 0.2089 -0.001948 NA
## Test set -12192.3 18058 15355 -1.513 1.933 0.3840 -0.632429 0.2873
mesesAPrever <- 6
Predição:
fcast_arima_6mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_6mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest6mth, col="red")
Acurácia:
accuracy(fcast_arima_6mth, tsMensalTest6mth)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 716.4 11391 8353 0.115 1.625 0.2089 -0.001948 NA
## Test set -10206.2 15249 12315 -1.258 1.539 0.3080 -0.585487 0.2658
mesesAPrever <- 12
Predição:
fcast_arima_1yr <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_1yr, include=48)
lines(tsMensalTest1yr, col="red")
Acurácia:
accuracy(fcast_arima_1yr, tsMensalTest1yr)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 716.4 11391 8353 0.1150 1.625 0.2089 -0.001948 NA
## Test set 4503.2 18305 15764 0.4428 1.841 0.3943 0.403569 0.4769
Técnica recente, publicada em 2011 (De Livera, Hyndman, and Snyder 2011).
A sigla TBATS é composta das iniciais das técnicas utilizadas: Trigonométrica, transformação Box-cox, correção de erro pela técnica ARMA, e componentes Trend e Sazonal.
Desenvolvido para séries temporais com características especiais de sazonalidade.
tbats_model <- tbats(tsDiarioTrain)
#tbats_model <- tbats(window(tsDiarioTrain, start=2005))
plot(tbats_model)
O método TBATS é indicado para a previsão do fluxo diário, pois:
Sazonalidades múltiplas: semanal e anual
Sazonalidade semanal é de alta frequência
diasAPrever <- 30
Predição:
tbats_fc30Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc30Days, include = 120)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(tbats_fc30Days, tsDiarioTest30Days)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set -13.13 742.6 436.1 -0.1878 2.095 0.2186 -0.0007912 NA
## Test set 412.21 689.8 601.6 1.3682 2.122 0.3015 0.6355962 0.3331
diasAPrever <- 45
Predição:
tbats_fc45Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc45Days, include = 120)
lines(tsDiarioTest45Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(tbats_fc45Days, tsDiarioTest45Days)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set -13.13 742.6 436.1 -0.1878 2.095 0.2186 -0.0007912 NA
## Test set 572.57 869.5 732.4 1.8836 2.518 0.3670 0.6215786 0.3952
Ver (R. J. Hyndman and Athanasopoulos 2013). O autor possui um livro exclusivamente sobre o assunto em (R. Hyndman et al. 2008).
Para cada um dos 15 métodos, há dois modelos possíveis: com erro aditivo e com erro multiplicativo.
ets(ts) sem argumentos além da série temporal ts determina automaticamente o método apropriado.ets apresenta o modelo escolhido e o AIC resultante.ets não trabalha com séries cuja sazonalidade é superior a 24 unidades temporais. Portanto, será utilizada somente pra prospecções realizadas com dados mensais (sazonalidade = 12).ets parametrizado automáticamente deve dar resultados muito precisos para prospecções de poucos pontos (por exemplo, 4 pontos).mesesAPrever <- 6
Modelo:
etsMensal <- ets(tsMensalTrain)
Prospecção:
fcastMensal6mth <- forecast(etsMensal, h=mesesAPrever)
Comparar predição com valores reais do conjunto de teste:
plot(fcastMensal6mth, xlab="Anos", ylab="Fluxo Mensal", include = 36);
lines(tsMensalTest6mth, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)
Acurácia:
accuracy(fcastMensal6mth, tsMensalTest6mth)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set 856.7 9849 7308 0.13394 1.442 0.1828 -0.02704 NA
## Test set -755.7 10937 8653 -0.06425 1.104 0.2164 -0.63165 0.2848
Há um fator multiplicativo no fator sazonal. Uma transformação Box-Cox demonstrou ser apropriada.
plot(tsMensalTrain, xlab="Anos", ylab="Fluxo Mensal")
Calcula-se o fator \(\lambda\) para a transformação:
lam <- BoxCox.lambda(tsMensalTrain)
Série Resultante:
tsMensalBoxCox <- BoxCox(x = tsMensalTrain, lam)
plot(tsMensalBoxCox, col="red")
mesesAPrever <- 24
Dados com sazonalidade removida:
tsMensal.stl <- stl(tsMensalTrain[,1], s.window=12)
# Seasonally adjusted data constructed by removing the seasonal component.
plot(seasadj(tsMensal.stl))
stlf_model <- stlf(tsMensalTrain[,1], lambda = lam)
stlf_fcast <- forecast(stlf_model, h=mesesAPrever)
plot(stlf_fcast, xlab="Anos", ylab="Fluxo Mensal")
lines(tsMensalTest2yr, col="red")
Acurácia:
accuracy(stlf_fcast, tsMensalTest1yr)
## ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
## Training set -25.42 8186 6175 -0.01584 1.238 0.1544 -0.02348 NA
## Test set 8132.24 17138 13696 0.92760 1.619 0.3425 -0.11417 0.4831
De Livera, Alysha M, Rob J Hyndman, and Ralph D Snyder. 2011. “Forecasting Time Series with Complex Seasonal Patterns Using Exponential Smoothing.” Journal of the American Statistical Association 106 (496): 1513–1527.
Hyndman, Rob J. 2014. “Forecasting Time Series Using R.” http://robjhyndman.com/talks/MelbourneRUG.pdf.
Hyndman, Rob J, and George Athanasopoulos. 2013. Forecasting: principles and Practice. OTexts.
Hyndman, Rob, Anne B. Koehler, J. Keith Ord, and Ralph D. Snyder. 2008. Forecasting with Exponential Smoothing: The State Space Approach (Springer Series in Statistics). 2008th ed. Springer. http://amazon.com/o/ASIN/3540719164/.
Leek, Jeffrey. 2014. “Practical Machine Learning - Lecture 27 - Forecasting.” Johns Hopkins Bloomberg School of Public Health. https://d396qusza40orc.cloudfront.net/predmachlearn/027forecasting.pdf.